// 2D SPDE model of thin film deposition model. See Ni Dong's dissertation for more details.
#include <cmath>
#include <fstream>
#include <iostream>
#include <string>
#include <iomanip>
using namespace std;
#include "model_SPDE2D_sincos.h"

// The following constant is used in update()
const double k_w = 3.3829e-12;
const double a_w = 0.6042;
const double E_w = 2.7e-3; // Unit = EV
const double k_c = 1.0274e-13;
const double a_c = 0.1669;
const double E_c = 1.9e-3; // Unit = EV
const double a_v = 15.55493;
const double k_v = 20.64504;
const double a_t = 0.02332;
const double k_t = 0.0261;
const double k_B = 8.617343e-5; // Boltzmann's constant, Unit = ev/K
// ---------------------------------------------------
model_SPDE2D_sincos::model_SPDE2D_sincos(int iMode,int iLatticeSize,double idt)
{
	int i_dim[2] = {1,1};
	int o_dim[4] = {(iMode+1)*(iMode+1),(iMode+1)*(iMode+1),(iMode+1)*(iMode+1),(iMode+1)*(iMode+1)};
	block::initialize(2,i_dim,4,o_dim,idt);

	int i,m,n;

	mode = iMode;
	LatticeSize = iLatticeSize;

	cov_z1 = block::outputs_next[0];
	cov_z2 = block::outputs_next[1];
	cov_z3 = block::outputs_next[2];
	cov_z4 = block::outputs_next[3];

	lambda = new double*[mode+1];
	for(m=0;m<mode+1;m++){
		lambda[m] = new double[mode+1];
	}

	for(i=0;i<4;i++){
		Km2[i] = new double[(mode+1)*(mode+1)];
	}	

	int index;
	for(m=0;m<=mode;m++){
		for(n=0;n<=mode;n++){
			index = index2D(m,n);
			if(m==0||2*m==LatticeSize||n==0||2*n==LatticeSize){
				Km2[0][index] = 0.0;
			}else{
				Km2[0][index] = 4.0/(pi*pi)*sin(pi*m/LatticeSize)*sin(pi*m/LatticeSize);
			}

			if(m==0){
				Km2[1][index] = 0.0;
			}else{
				Km2[1][index] = 4.0/(pi*pi)*sin(pi*m/LatticeSize)*sin(pi*m/LatticeSize);
			}

			if(m==0||2*m==LatticeSize){
				Km2[2][index] = 0.0;
			}else{
				Km2[2][index] = 4.0/(pi*pi)*sin(pi*m/LatticeSize)*sin(pi*m/LatticeSize);
			}

			if(n==0 || 2*n==LatticeSize || m==0){
				Km2[3][index] = 0.0;
			}else{
				Km2[3][index] = 4.0/(pi*pi)*sin(pi*m/LatticeSize)*sin(pi*m/LatticeSize);
			}
		}
	}

	reset();
}
// ---------------------------------------------------
model_SPDE2D_sincos::~model_SPDE2D_sincos()
{
	int m;
	for(int i=0;i<4;i++)
	{
		delete []Km2[i];
	}

	for(m=0;m<mode+1;m++){
		delete []lambda[m];
	}
	delete []lambda;
}
// ---------------------------------------------------
void model_SPDE2D_sincos::reset(){
	block::time = dt;
	int m,n;
	int index;

	mean_z200 = 0.0;
	for(m=0;m<mode+1;m++){
		for(n=0;n<mode+1;n++){
			index = index2D(m,n);
			cov_z1[index] = 0.0;
			cov_z2[index] = 0.0;
			cov_z3[index] = 0.0;
			cov_z4[index] = 0.0;
			block::outputs[0][index] = 0.0;
			block::outputs[1][index] = 0.0;
			block::outputs[2][index] = 0.0;
			block::outputs[3][index] = 0.0;
			block::outputs_next[0][index] = 0.0;
			block::outputs_next[1][index] = 0.0;
			block::outputs_next[2][index] = 0.0;
			block::outputs_next[3][index] = 0.0;
		}
	}

	T = 610;
	W = 0.5;
	c = eq_c(T,W);
	c2 =eq_c2(T,W);
	sigma = sqrt(eq_sigma2(T,W));
	//c = 0.5;
	//c2 = 0.1;
	//sigma = 0.1;

	for(m=0;m<mode+1;m++){
		for(n=0;n<mode+1;n++){
			lambda[m][n] = -4*c2*(m*m+n*n);
		}
	}
}
// ---------------------------------------------------
void model_SPDE2D_sincos::update(double sysTime, double iT, double iW,const double *z_mean0){
	if(time>sysTime+dt) return;

	int m,n;
	int index;
	double c_200;
	double temp1;
	double temp2;
	if(time<=sysTime+dt){
		block::output();
	}
	while(time<=sysTime+dt){
		if(iT!=T || iW!=W){
			T = iT;
			W = iW;
			c = eq_c(T,W);
			c2 = eq_c2(T,W);
			sigma = sqrt(eq_sigma2(T,W));
			for(m=0;m<mode+1;m++){
				for(n=0;n<mode+1;n++){
					lambda[m][n] = -4*c2*(m*m+n*n);
				}
			}
		}
//		c  = 0.5;
//		c2 = 0.1;
//		sigma = 0.1;
		c_200 = LatticeSize*c;
		mean_z200 = mean_z200+c_200*dt;
		for(m=0;m<=mode;m++){
			for(n=0;n<=mode;n++){
				index = index2D(m,n);
				temp1 = exp(2.0*lambda[m][n]*dt);
				temp2 = sigma*sigma*(temp1-1.0)/(2.0*lambda[m][n]);
				/*if(m==0&& n==0){
					cov_z2[index] = cov_z2[index]+sigma*sigma*dt;
				}

				if(m!=0&&n!=0){
					cov_z1[index] = cov_z1[index]*temp1+temp2;
					cov_z2[index] = cov_z2[index]*temp1+temp2;
					cov_z3[index] = cov_z3[index]*temp1+temp2;
					cov_z4[index] = cov_z4[index]*temp1+temp2;
				}

				if(m==0 && n!=0){
					cov_z2[index] = cov_z2[index]*temp1+temp2;
					cov_z4[index] = cov_z4[index]*temp1+temp2;
				}

				if(m!=0 && n==0){
					cov_z2[index] = cov_z2[index]*temp1+temp2;
					cov_z3[index] = cov_z3[index]*temp1+temp2;
				}*/

				if(m>0&&n>0){
					cov_z1[index] = cov_z1[index]*temp1+temp2;
				}

				if(m==0&&n==0){
					cov_z2[index] = cov_z2[index]+sigma*sigma*dt;
				}else{
					cov_z2[index] = cov_z2[index]*temp1+temp2;
				}

				if(m>0){
					cov_z3[index] = cov_z3[index]*temp1+temp2;
				}

				if(n>0){
					cov_z4[index] = cov_z4[index]*temp1+temp2;
				}
			}
		}
		time += dt;
	}
	return;
}
// ---------------------------------------------------
void model_SPDE2D_sincos::update(double sysTime)
{
	double T = inputs[0][0];
	double W = inputs[1][0];
	update(sysTime,T,W);
}
// ---------------------------------------------------
double model_SPDE2D_sincos::getR2()
{
	int m,n;
	int index;
	double r2 = 0.0;

	for(m=0;m<=mode;m++){
		for(n=0;n<=mode;n++){
			if(m!=0||n!=0){
				index = index2D(m,n);
				r2+=outputs[0][index]+outputs[1][index]+outputs[2][index]+outputs[3][index];
				//r2+= cov_z2[index];
			}
		}
	}

	//r2 = cov_z2[index2D(0,1)];
	r2 /= (pi*pi); // eq(5.34)
	return r2;
}
// ---------------------------------------------------
double model_SPDE2D_sincos::getM2()
{
	int m,n;
	int index;
	double m2 = 0.0;
	for(m=0;m<=mode;m++){
		for(n=0;n<=mode;n++){
			index = index2D(m,n);
			m2+= Km2[0][index]*outputs[0][index]
			+Km2[1][index]*outputs[1][index]
			+Km2[2][index]*outputs[2][index]
			+Km2[3][index]*outputs[3][index];
		}
	}
	return m2;
}
// ---------------------------------------------------
double model_SPDE2D_sincos::getH()
{
	return sqrt(cov_z2[index2D(0,0)])/LatticeSize;
}
// ---------------------------------------------------
double model_SPDE2D_sincos::eq_c(double T,double W)
{
	return W*(1-k_w*pow(W,-a_w)*exp(k_B*T/E_w));
}

double model_SPDE2D_sincos::eq_c2(double T,double W)
{
	return k_c*exp(k_B*T/E_c)*pow(W,-a_c)/(LatticeSize*LatticeSize);
}

double model_SPDE2D_sincos::eq_sigma2(double T,double W)
{
	double index = (a_t+k_t*W)*T-a_v-k_v*W;
	return pi*pi/(LatticeSize*LatticeSize)*(1+exp(index));
}
